##Table 3: COVID EIDL##

EIDL<-read.csv("EIDL_businessonly.csv");

##modify the setting to add additional parameters in the estimation of counterfactual distribution##

trace(prep_data_for_fit, edit=TRUE);

##replace  

"      if (sum(!is.na(rn)) > 0) {
        rn <- rn[order(rn)]
        rn_vector <- c()
        for (i in rn) {
            rn_dummy_name <- paste0("rn_", i)
            data_binned[[rn_dummy_name]] <- ifelse(data_binned$bin%%i == 
                0, 1, 0)
            rn_vector <- c(rn_vector, rn_dummy_name)
        }
        if (length(rn) == 2 & (rn[2]%%rn[1] == 0)) {
            data_binned[[rn_vector[1]]] <- ifelse(data_binned[[rn_vector[1]]] == 
                1 & data_binned[[rn_vector[2]]] == 1, 0, data_binned[[rn_vector[1]]])
        }
    }
    else {
        rn_vector <- ""
    }
    rhs_vars <- c("zstar", extra_fe_vector, polynomial_vector, 
        rn_vector, bins_excluded_all)" by
"
 if (sum(!is.na(rn)) > 0) {
        rn <- rn[order(rn)]
        rn_vector <- c()
        for (i in rn) {
            rn_dummy_name <- paste0("rn_", i)
            data_binned[[rn_dummy_name]] <- ifelse(data_binned$bin%%i == 
                0, 1, 0)
            rn_vector <- c(rn_vector, rn_dummy_name)
        }
        rn_dummy_name1 <- paste0("before_5000rn")
        data_binned[[rn_dummy_name1]] <- ifelse((data_binned$bin + 
            1000)%%5000 == 0, 1, 0)
        rn_dummy_name2 <- paste0("before_10000rn")
        data_binned[[rn_dummy_name2]] <- ifelse((data_binned$bin + 
            1000)%%10000 == 0, 1, 0)
        rn_vector <- c(rn_vector, rn_dummy_name1, rn_dummy_name2)
        data1 <- cbind(data_binned[[rn_vector[1]]], data_binned[[rn_vector[2]]], 
            data_binned[[rn_vector[3]]], data_binned[[rn_vector[4]]])
        data2 = as.data.frame(t(apply(data1, 1, function(x) {
            x = ifelse(duplicated(x, fromLast = T), 0, x)
        })))
        data_binned[[rn_vector[1]]] <- data2[, 1]
        data_binned[[rn_vector[2]]] <- data2[, 2]
        data_binned[[rn_vector[3]]] <- data2[, 3]
        data_binned[[rn_vector[4]]] <- data2[, 4]
    }
    if (sum(!is.na(rn)) > 0) {
        rn <- rn[order(rn)]
        rninter_vector <- c()
        for (i in 1:(length(rn) + 2)) {
            rninter_dummy_name <- paste0("rninter_", i)
            data_binned[[rninter_dummy_name]] <- data_binned[[rn_vector[i]]] * 
                data_binned$z_rel
            rninter_vector <- c(rninter_vector, rninter_dummy_name)
        }
    }
    rhs_vars <- c("zstar", extra_fe_vector, polynomial_vector, 
        rn_vector, rninter_vector, bins_excluded_all)
"
##

a=24000;
b=40500;
c=0;
zstar=25000;
x0=500;
binwidth=x0;
y0=5;


binned_data <- bunching::bin_data(z_vector = EIDL$facevalueofdirectloanorloanguara,zstar =zstar, binwidth = x0,
                          bins_l =a/x0, bins_r = b/x0);


bins_excl_r <- 1
         
            firstpass_prep <- bunching::prep_data_for_fit(binned_data,zstar = zstar, binwidth = x0,
                          bins_l =a/x0, bins_r = b/x0, rn=c(5000,10000),bins_excl_l=c/x0,bins_excl_r=bins_excl_r,poly=y0)
            firstpass <- bunching::fit_bunching(firstpass_prep$data_binned, 
                firstpass_prep$model_formula,binwidth=x0, notch=TRUE)
            B_below <- firstpass$B_zl_zstar
            M_above <- -firstpass$B_zstar_zu
                       if (M_above > B_below) {
                stop("Missing mass above zstar is larger than bunching mass below. Are you sure this is a notch?")
            }
            available_bins <- (max(firstpass_prep$data_binned$bin) - 
                zstar)/binwidth
            DoF_remaining <- nrow(firstpass_prep$data_binned) - 
                nrow(firstpass$coefficients) - 1
            notch_iterations_bound <- min(available_bins, DoF_remaining)
            zu_bin <- bins_excl_r
            tmp_firstpass_prep <- firstpass_prep
            while ((B_below > M_above) & (zu_bin < notch_iterations_bound)) {
                zu_bin_final <- zu_bin
                zu_bin <- zu_bin + 1
                newvar <- paste0("bin_excl_r_", zu_bin)
                tmp_firstpass_prep$data_binned[[newvar]] <- ifelse(tmp_firstpass_prep$data_binned$z_rel == 
                  zu_bin, 1, 0)
      tmp_firstpass_prep$model_formula <- stats::as.formula(paste(Reduce(paste, 
                  deparse(tmp_firstpass_prep$model_formula)), 
                  newvar, sep = " + "))

   tmp_firstpass <- bunching::fit_bunching(tmp_firstpass_prep$data_binned, 
                  tmp_firstpass_prep$model_formula,binwidth=x0, notch=TRUE)

                B_below <- tmp_firstpass$B_zl_zstar
                M_above <- -tmp_firstpass$B_zstar_zu
                if (B_below > M_above) {
                  firstpass_prep <- tmp_firstpass_prep
                  firstpass <- tmp_firstpass
                  zu_bin_final <- zu_bin
                }
            }
     
mb<-zu_bin*x0+25000

  B<- tmp_firstpass$B_zl_zstar/sum(binned_data$freq)
      
 forbootpass1<- tmp_firstpass;
forbootprep1<- tmp_firstpass_prep;

 ratio=(1-zstar/(mb));
beta<-1/3;
alpha1<-1/3;
R<-1+0.0375
    cost1=((1/alpha1)*(1-(1-ratio)^(alpha1))-ratio)*R;
nu=0.83
alpha2<-beta*nu/(1-(1-beta)*nu)
cost2=((1/alpha2)*(1-(1-ratio)^(alpha2))-ratio)*R;


n_boot = 1000

    # retrieve data and model from firstpass_prep
     data_for_boot <- forbootprep1$data_binned
   residuals <-forbootpass1$residuals;

    # get vector of bootstrapped betas
    boot_results <- lapply(seq(1:n_boot), function(i) {
        # adjust frequencies using residual
        data_for_boot$freq <- data_for_boot$freq_orig+  sample(residuals, replace = TRUE)
        # make this "freq" so we can pass to fit_bunching which requires "freq ~ ..."
      
 bins_excl_r <- 1

  boot_firstpass_prep <- bunching::prep_data_for_fit(data_for_boot,zstar = zstar, binwidth = x0,
                          bins_l =a/x0, bins_r = b/x0, rn=c(5000,10000),bins_excl_l=c/x0,bins_excl_r=bins_excl_r,poly=y0)
            boot_firstpass <- bunching::fit_bunching(boot_firstpass_prep$data_binned, 
                boot_firstpass_prep$model_formula,binwidth=x0, notch=TRUE)
    
    # extract bunching mass below and missing mass above zstar
    B_below <-   boot_firstpass$B_zl_zstar
    M_above <- -boot_firstpass$B_zstar_zu
 
    # if(M_above > B_below) == TRUE, we start shifting zu bins up until B = M.
    # first, must set upper bound on number of iterations.
    # bins above zstar cannot be more than min of:
    #   1. available bins
    #   2. remaining DoF

       if (M_above > B_below) {
  Bestimate<-   B_below;
    booted_zU_notch=0;
    ratio=(1-zstar/(zstar+booted_zU_notch*x0));
     cost1=((1/alpha1)*(1-(1-ratio)^(alpha1))-ratio)*R;
cost2=((1/alpha2)*(1-(1-ratio)^(alpha2))-ratio)*R;


tempresult<-cbind(booted_zU_notch,Bestimate,ratio,cost1,cost2)
     } else{                      
    available_bins <- (max(boot_firstpass_prep$data_binned$bin) - zstar)/binwidth
    DoF_remaining <- nrow(boot_firstpass_prep$data_binned) -  nrow(boot_firstpass$coefficients) - 1
    notch_iterations_bound <- min(available_bins, DoF_remaining)
    
    # start with first bin above being bins_excl_r = 1
    zu_bin <- bins_excl_r
    
    # create temporary object for firstpass_prep (we will be editing this as we expand the window for zu_bin)
    tmp_firstpass_prep <- boot_firstpass_prep
    
    while ((B_below > M_above) & (zu_bin < notch_iterations_bound)) {
        # keep previous zu_bin in case next gives us B > M
        zu_bin_final <- zu_bin
        # now try expanding by 1 bin
        zu_bin <- zu_bin + 1
        # add next bin_excl_r dummy as column to data
        newvar <- paste0("bin_excl_r_", zu_bin)
        tmp_firstpass_prep$data_binned[[newvar]]  <- ifelse(tmp_firstpass_prep$data_binned$z_rel == zu_bin,1,0)
        # add next order bin_excl_r to formula
         tmp_firstpass_prep$model_formula <- stats::as.formula(paste(Reduce(paste, deparse(tmp_firstpass_prep$model_formula)), newvar, sep = " + "))
        # re-fit model using the now expanded zu
        tmp_firstpass <- bunching::fit_bunching(tmp_firstpass_prep$data_binned, tmp_firstpass_prep$model_formula,binwidth=x0,  notch=TRUE)
        # get new B below and M above
        B_below <- tmp_firstpass$B_zl_zstar
        M_above <- -tmp_firstpass$B_zstar_zu
        
        # if M is not larger, update firstpass with this (last iterations of loop will have M > B)
        if(B_below > M_above) {
            boot_firstpass_prep <- tmp_firstpass_prep
            boot_firstpass <- tmp_firstpass
            zu_bin_final <- zu_bin
        }
        
    }

  firstpass1<- tmp_firstpass;
firstpass_prep1<- tmp_firstpass_prep;

    # assign final zu_bin_final to bins_excl_r (used for plotting)


      

   Bestimate<- firstpass1$B_zl_zstar;
Mestimate<-firstpass1$B_zstar_zu*(-1);
    booted_zU_notch=zu_bin;

    ratio=(1-zstar/(zstar+booted_zU_notch*x0));
     cost1=((1/alpha1)*(1-(1-ratio)^(alpha1))-ratio)*R;
cost2=((1/alpha2)*(1-(1-ratio)^(alpha2))-ratio)*R;


tempresult<-cbind(booted_zU_notch,Bestimate,ratio,cost1,cost2)
}

   
        return(tempresult)
    })



    zU_Notch_boot <-do.call(rbind,boot_results)
zU_Notch_boot<-as.data.frame(zU_Notch_boot);
zU_Notch_boot<-subset(zU_Notch_boot,zU_Notch_boot$dif>=0);

zU_Notch_boot<-subset(zU_Notch_boot,zU_Notch_boot$booted_zU_notch<=50);



zU_bin_sd<-sd(zU_Notch_boot$booted_zU_notch);
B_sd<-sd(zU_Notch_boot$Bestimate)/sum(binned_data$freq);
mb_sd<-zU_bin_sd*x0
ratio_sd<-sd(zU_Notch_boot$ratio);
cost1_sd<-sd(zU_Notch_boot$cost1);
cost2_sd<-sd(zU_Notch_boot$cost2);


result<-cbind(B,B_sd,mb,mb_sd,ratio,ratio_sd,cost1,cost1_sd,cost2,cost2_sd);

result


